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Abstract —Stochastic economic dispatch models address uncer¬ 
tainties in forecasts of renewable generation output by consid¬ 
ering a finite number of realizations drawn from a stochastic 
process model, typically via Monte Carlo sampling. Accurate 
evaluations of expectations or higher-order moments for quanti¬ 
ties of interest, e.g., generating cost, can require a prohibitively 
large number of samples. We propose an alternative to Monte 
Carlo sampling based on Polynomial Chaos expansions. These 
representations are based on sparse quadrature methods, and 
enable accurate propagation of uncertainties in model param¬ 
eters. We also investigate a method based on Karhunen-Loeve 
expansions that enables us to efficiently represent uncertainties in 
renewable energy generation. Considering expected production 
cost, we demonstrate that the proposed approach can yield 
several orders of magnitude reduction in computational cost 
for solving stochastic economic dispatch relative to Monte Carlo 
sampling, for a given target error threshold. 

Index Terms —Stochastic Economic Dispatch, Monte Carlo 
Sampling, Polynomial Chaos Expansion, Karhunen-Loeve Ex¬ 
pansion. 

I. Introduction 

U NIT commitment (UC) is the fundamental process of 
scheduling thermal generating units in advance of op¬ 
erations in the electric power grid [1], The objective is to 
minimize overall production costs to satisfy forecasted load 
for electricity, while respecting operational constraints of both 
transmission elements (e.g., thermal limits) and generators 
(e.g., ramping limits). Economic dispatch (ED) is a closely 
related operations problem, in which cost minimization is 
performed to identify an optimal set of power output levels 
for a fixed set of on-line thermal generating units. Typically, 
UC and ED are respectively formulated as mixed-integer 
and linear optimization problems, and routinely solved using 
commercial solvers. Despite recent improvements in load 
forecasting technology, next-day predictions are imperfect, 
with errors on average in the 1-3% range and exceeding 10% 
on specific days [2], To account for such inaccuracies, reserve 
margins are universally imposed in operations planning. These 
margins implicitly deal with uncertainty in load forecasts, by 
ensuring there is sufficient generation capacity available to 
meet unexpectedly high load during operations. 

An alternative approach to dealing with forecast errors is 
to explicitly model the uncertainty, typically via a finite set 
of sampled realizations from a stochastic process model. This 
approach results in a stochastic ED model (SED), in which 


the objective typically is to minimize the expected production 
cost across a set of load scenarios [3], [4]. By explicitly 
representing the inherent uncertainty in load forecasts, a 
SUC solution ensures sufficient flexibility to meet a range of 
potential realizations during next-day operations. Further, by 
explicitly representing uncertainty, reliance on reserve margins 
is reduced, yielding less costly solutions than those obtained 
under deterministic ED models. Increasing penetration levels 
of renewables (e.g., wind and solar) generation accentuate the 
differences between stochastic and deterministic grid opera¬ 
tion problems, due to increased errors in next-day forecasts; 
relative to load, accurate prediction of next-day meteorological 
conditions is very difficult. 

Despite the conceptual appeal of stochastic variants of grid 
operations models, they are not yet used in in practice due to 
their well-known computational difficulty [5]. This difficulty is 
primarily driven by the number of forecast samples required 
to achieve high-quality, robust solutions to models such as 
stochastic UC and ED. 

Uncertainties such as those found in stochastic UC are ubiq¬ 
uitous in both power systems operations and planning, and the 
importance of credibly accounting for them is well-recognized. 
In their seminal work, Takriti et al. introduced a model and 
solution technique for the problem of generating electric power 
when loads are uncertain [4]. More recent work on addressing 
grid operations under uncertainty has focused on various 
sources of uncertainty and different frameworks for modeling 
uncertainty (e.g., stochastic programming, chance constraints, 
and robust optimization). For example, [6] considers a stochas¬ 
tic UC in which the availability of consumer demand response 
(DR) is uncertain. DR uncertainty was modeled using a set of 
scenarios, and a chance constraint was imposed to ensure loss- 
of-load probability is below a pre-defined risk level. Chen el al. 
[7] considered the combined UC and ED problem under both 
random and targeted component failures, where allowable loss- 
of-load was parameterized by the contingency size. Bertsimas 
et al. [8] proposed a two-stage adaptive robust UC model given 
uncertainty in nodal net injections, and developed a solution 
approach based on a combination of Benders decomposition 
and outer approximation. 

Despite recent advances, the general lack of advanced 
methods to accurately model and represent uncertainty and the 
inability of scenario-based approaches to solve industrial sized 
stochastic optimization problems have led researchers to seek 
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alternatives, in both grid operations and planning contexts. 
For example, Thiam and DeMarco [9] argue: “ Simply put, 
when uncertainty is credibly accounted for such methods yield 
solutions for economic benefit of a transmission expansion 
in which the “error bars” are often larger than the nominal 
predicted benefit In other words, only a limited number of 
samples could be considered while sustaining computational 
tractability, which in turn impacts the ultimate utility of a 
solution. Of course, it is not possible to change the nature of 
uncertainties, such that if uncertainties are so large that they 
fail to provide useful information. However, as we demonstrate 
in our computational experiments, it is possible to significantly 
reduce the impact of errors introduced by modeling and 
sampling of uncertainty. 

In this paper, we adopt advanced modeling and sampling 
techniques from the uncertainty quantification (UQ) commu¬ 
nity, and leverage them to impact power systems operations 
problems such as stochastic UC and ED. Such techniques 
have been successfully applied in many areas of computational 
science and engineering, with significant success (see e.g., [10] 
and references therein). The need for accurate estimation 
of uncertain model outputs, along with the prohibitive cost 
of requisite large numbers of Monte Carlo (MC) samples, 
have led to the development of more efficient alternatives. 
Specifically, we employ Polynomial Chaos expansions [11] 
to represent uncertain model inputs in terms of sets of or¬ 
thogonal polynomials of standard random variables. The task 
of propagating this functional representation to model outputs 
can be achieved via several pathways. In our experiments, 
we employ a Galerkin projection technique in conjunction 
with sparse quadrature methods. We demonstrate that our 
approach yields a one to two order of magnitude reduction 
in the number of samples (scenarios) required to estimate 
expected production cost.relative to MC, depending on given 
target error thresholds. Consequently, our approach has the 
potential to dramatically reduce the computational difficulty 
of stochastic grid operations problems, significantly reducing 
a major barrier to their use in practice. 

The remainder of this paper is organized as follows. We 
briefly introduce our stochastic ED formulation in Section II 
to provide context for our research. In Section III we describe 
a method to efficiently model the uncertainties in renewables 
power production, with an emphasis on wind generation. 
Section IV reviews key concepts in the representation of 
uncertainty using Polynomial Chaos and discusses the use 
of our surrogate models of renewable resource output for 
stochastic ED. We then empirically analyze the accuracy of 
our surrogates on standard IEEE test problems in Section IV-C, 
and conclude in Section V. 

II. Stochastic Economic Dispatch 

E abstractly denote the set of unit commitment con¬ 
straints (i.e., operational and physical constraints on 
physical units) as X and let x £ X denote a vector of (binary) 
unit commitment decisions. In the context of this paper, we 
assume that the unit commitment decisions x are fixed; thus, 
x are parameters in our SED model. 


We treat renewable power generation as a random field and 
present an efficient approach to represent these random fields 
in subsequent sections. In principle, one can also consider 
uncertain loads and similarly represent them as random fields 
[12], For conciseness, we focus strictly in our experiments 
on uncertain renewable outputs and without loss of generality 
consider deterministic loads. 

Because uncertain renewables generation is represented 
using random fields and ultimately, as will be illustrated 
below, as functions of a vector of random variables £(co), 
the corresponding production cost 2(jr,£({D)) is similarly 
uncertain or random. The expected production cost, denoted 
Q(x), is defined as 

Q(x)=EsQ(x£(< 0 )) ( 1 ) 

and the (multi-period) stochastic economic dispatch problem 
under a fixed unit commitment x is given by 

fi(*>$(«>)) =, “ in >0 e £ £ +E E Mc li ( 2 ) 

f,p> 0,9>0,e teTgeG teTieN 

s.t. 

Y J P , Mu))+ £ p' g + £ f e -£ vu 

r£Rj g&Gj e ^E j e^Ei 

(3a) 

Be(%-%)-£= 0, Ve = (i,j),t (3b) 

F e <f e <F e , Vej (3c) 

Pgx'g <Pg< Pgx'g, Vg, t (3d) 

p'g-p^ 1 <R u g xf g l +S u g {xf g -tf~ l )+P g (\-x‘ g ), Vg,t (3e) 

p'g 1 -p' g <R d g ^ g +S d g {f~ l -x t g )+P g (l-tf g ~ l ), Vgf (3D 

where R" (R d g ) and S g (S'!,) represent nominal ramp-up (ramp- 
down) and startup (shutdown) rates, respectively. 

Here, subscripts i and j denote bus indexes defined over 
bus set N, while the superscript t denotes specific time 
periods in the planning horizon t = Subscript g 

denotes the generator index defined over generator set S, while 
the subscript e = ( i,j ) denotes the line index (and terminal 
buses (i,j)) defined over transmission line set £. Renewable 
generation power (a parameter) is denoted by p' r while power 
from thermal generators (a variable) is denoted by p' The 

decision variables q\ denotes the load shedding quantity at bus 
i at time period t. Because renewable generation is a function 
of RVs, solution variables are necessarily RVs. The output of 
interest Q is a function of the same set of random variables 
that are used to define the uncertain renewables. 

The optimization objective in SED is to minimize the 
expected total production and loss-of-load costs. The first 
term in Eq. (2) represents total production cost, while the 
second term represents the loss-of-load penalty. The load¬ 
shedding penalty employs a large positive number M, typically 
$5000 or $6000 [5]. Eq. (3) specifies operational and physical 
constraints on grid components based on a direct current (DC) 
power flow model, and includes (in order): power balance at 
each time period and bus (3a); the power flow f e through 
each line e = ( i,j ) as a function of voltage phase angle 
differences of the terminal buses, (0, —0 ; ), and the reciprocical 
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of the line reactance, B e (3b); lower and upper bounds 
(F_g and F e , respectively) on line power flow (3c); lower 
and upper bounds (P„ and P g , respectively) for the standard 
generator power (3d); and generation ramp-up and ramp-down 
constraints for pairs of consecutive time periods (3e, 3f). 

A quadratic production cost function is often employed in 
economic dispatch models, e.g., as follows: 

c g(P t g) =a g*?g+ b gP t g + c g(P t g ) 2 - ( 4 ) 

Equation (4) can be accurately approximated using a set of 
piecewise linear segments. For conciseness, we omit these 
standard linearization steps. For details concerning lineariza¬ 
tion of the quadratic production cost functions, see [1], 
Within the broader context of combined stochastic UC&ED, 
the stochastic ED model is embedded as a sub-problem in 
the UC model. In a combined stochastic UC and ED model, 
the first-stage decisions are the unit commitment selections 
x , while the optimization objective is to minimize expected 
production costs. In the second (recourse) decision stage, 
uncertain renewable production results in uncertain (sample- 
dependent) recourse decisions for the dispatch and load¬ 
shedding variables, and consequently uncertain production and 
load-shedding costs. First-stage unit commitment decisions are 
determined by taking their future impacts into consideration. 
These future impacts are quantified by the recourse function 
Q{x), which computes the expected value of production cost 
for a given (fixed) unit commitment x. 

Next, we describe a method for representing uncertain wind 
generation as a spectral expansion that decouples the deter¬ 
ministic space (time) from the stochastic space and provides 
an avenue to significantly reduce the dimensionality of the 
stochastic space. 

III. Modeling Wind Uncertainty via 
Karhunen-Loeve Expansions 
Accurate, efficient, and low-dimensional representations 
of uncertainties are essential for the success of stochastic 
grid operations models. We now discuss how to construct a 
lower dimensional representation of wind power uncertainty. 
Towards this goal, we explore representing uncertain wind 
production time profiles p' r via Karhunen-Loeve (KL) expan¬ 
sions [13]. We start in Section III-A by considering data from 
NREL’s Western Wind Dataset [14], to assess the feasibility 
of using KL expansions for representing uncertainty in wind 
speed profiles. We then propose, in Section III-B, a model 
for generating wind power samples that is consistent with 
uncertainties observed in current forecast models. We avoid 
representing uncertainty in wind power directly, due to the 
discrete nature of the cut-off threshold. 

A. Karhunen-Loeve Expansions 

The KL expansion represents stochastic processes by a 
linear combination of orthogonal modes. In order to assess 
the feasibility of this approach for representing wind gener¬ 
ator output, we consider data for three wind sites extracted 
from NREL’s Western Wind Dataset [14], Two of these 
sites are located in Wyoming and are geographically close: 



W [m/s] 


Fig. 1. Rated power output P vs. wind speed W at site #15414 for year 
2004. Blue dots represent instantaneous measurements while the red curve 
represents power data filtered with a top-hat filter of width Aw =0.05 [m/s], 
and interpolated via cubic splines. 


site #15414 located at (41.48/V, 105.14W) and site #16238 
located at (41.61 A'. 105.1 IW). The third site, #3560 located at 
(35.08A, 118.41VT) in California, was selected to be far from 
the first two sites. The proximity of the first two sites allows 
us to explore correlations in wind power, while the third site is 
sufficiently distant from the first two sites such that no spatial 
correlations in their production exist. 

For each site, wind speed and power data is available at 10 
minute intervals for the years 2004 through 2006. Fig. 1 shows 
a scatter plot of rated power output vs. wind speed at wind 
site #15414. The rated power output is zero for wind speeds 
smaller than a threshold value of approximately 3.2 m/s. At 
speeds greater then approximately 26 m/s the generators are 
turned off, for safety reasons. 

We take the following steps to post-process the NREL wind 
data into wind samples for use in our SED model: 

1) Construct hourly averages for the wind speed; consider 
wind speed data for each day to be an independent 
sample from a 24-dimensional random field. 

2) Choose a range of dates and assemble a set of 24- 
dimensional samples from this date range. For example, 
considering the month of January for 2004 through 2006 
leads to 93 samples. This step was adopted to account 
for seasonal changes in wind patterns. 

Fig. 2 shows select daily log-transformed (Wl = logfH')) wind 
speed samples at site #15414. We adopted this transformation 
to ensure positivity of wind speed samples, generated via the 
algorithm described below. 

We represent Wl via a KL expansion 

W L {t, to) = {W L (t, CO)) + £ Vhfk(t)^k (®) (5) 

k=l 

where (VL,(f,C0)) denotes the mean of Wj,(f,G)), fk(t) and 
Xl are respectively eigenfunctions (eigenvectors in a discrete 
setting) and eigenvalues of the covariance matrix Eyr,. of 
Wl(p to), and ^ denotes uncorrelated random variables with 
zero mean and unit variance. Projecting realizations of Wl 
onto fic leads to samples of These samples are generally 
not independent. In the special case where Wl is a Gaussian 
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Fig. 2. Select daily wind data from January 2004-2006 at site #15414. The 
red curve corresponds to the mean of these daily samples. 


Fig. 4. Percentages of the total variance explained by truncated KL 
expansions with an increasing number of modes, for our three wind sites. 



Fig. 3. The mean, and first three eigenvectors, /i, / 2 , and fi, at 

several wind sites based on data from January 2004-2006. 

random process, ^ are i.i.d. standard normal random variables. 

If known, the covariance matrix E w can be specified analyti¬ 
cally. Otherwise, if sufficient samples are available, the covari¬ 
ance matrix can also be estimated from these realizations. For 
In our experiments, we estimate the covariance matrix using 
daily samples for select times of the year during 2004-2006. 
Once the covariance matrix is available, the eigenvalues and 
eigenvectors in Eq. (5) are given by solutions of the Fredholm 
equation of second kind 

J Z WL (t,s)f(s)ds = Xf(t). 

We discretize and solve this equation using the Nystrom 
method [15], considering a mid-point quadrature rule to in¬ 
tegrate over the discrete 24-hour period. 

Fig. 3 shows the mean and first few KL modes correspond¬ 
ing to the month of January at the three wind sites considered 
in this study. The large degree of similarity for the KL modes 
between the three sites suggest a similar structure for the time 
correlations in the wind speed. We will further explore the 
structure of the covariance matrices E w in Section III-B. 

Because the stochastic dimensions are uncorrelated, the total 



Fig. 5. Reconstruction of daily wind samples via truncated KL expansion. 
Left plot corresponds to Jan 1, 2004 and the right plot to Jan 26, 2004. 


variance of the stochastic process is given as the sum of 
variances from individual KL mode, as follows: 

24 

Var[W L (f,C 0 )] = £ \ k fi(t)Var[Q . 

” 1 =1 

Given that eigenvectors are orthonormal with respect to the 
deterministic space (the discretized time axis in our study), it 
follows that: 

r 24 

/ Vai\W L (t,(o)\dt= ( 6 ) 

Jt k—\ 

As a result, an N-truncated expansion, N < 24, will explain 
100 x j 

of the total variance of the random field. Fig. 4 shows the 
dependence of the fractional variance given in Eq. (7) on the 
number of terms N in the truncated KL expansions at the three 
wind sites. It is evident that for all sites, N = 6 modes are 
sufficient to capture approximately 95% of the total variance 
in the daily wind profiles. We also examine the re-construction 
of select daily wind samples with truncated KL expansions, 
shown in Fig. 5. The results indicate that N = 6 KL modes 
are sufficient to represent most of the daily variability. Similar 
results are observed for other, randomly selected, samples. 

Next, we test the degree of dependence between KL random 
variables. We start by constructing empirical cumulative dis¬ 
tribution functions (CDFs) from RV samples. These samples 
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Fig. 6 . Empirical CDFs (depicted in gray) for through ^15 at (a) site 
16238 and (b) site 3560. The CDF for a standard normal random variable is 
shown in blue, for reference. 


TABLE I 

Distance correlation factors for £1 - i; 4 between wind sites. 


RV 

16238-15414 

16238-3560 

15414-3560 


0.91 

0.24 

0.24 


0.84 

0.13 

0.18 


0.67 

0.18 

0.19 


0.53 

0.23 

0.22 


are obtained by projecting each daily wind speed sample onto 
the corresponding KL mode. Fig. 6 shows the empirical CDFs 
for through ^15 at two select wind sites. Visual inspection 
of these results and comparison with the CDF of a standard 
normal RV indicate strong similarities. Because these RVs 
are uncorrelated by construction, modeling them as standard 
normals implies that they are also independent. 

We further investigate relationships among the standard 
normal RVs across the sites we consider in our experiments. 
To this end we employ distance correlation factors [16], [17] 
between pairs of RVs corresponding to the same KL mode 
at different sites. These factors are shown in Table I. Smaller 
values, towards zero, indicate negligible dependence between 
pairs of RVs, while larger values, close to 1, indicate a 
strong dependence. These results indicate a strong correlation 
between sites #15414 and #16238 for the first two KL modes, 
while the same KL mode at site #3560 shows little correlation 
with the first two sites. The third and fourth RVs, as well as 
the other RVs (not shown), show little correlation between 
sites. This is somewhat to be expected given the nature of 
turbulence. Specifically, low KL modes are associated with 
large scale structures which are likely similar if sites are 
geographically situated, while higher order KL modes are 
associated with smaller scale structures with much faster eddy 
turnover times. 

In experiments described below, we consider a stochastic 
space with N = 6 dimensions at each site. Given that we 
model the first two modes at two sites as dependent, the 
total dimensionality of the stochastic space is 3x6 — 2 = 16. 
Given the representation of daily wind profiles as truncated 
KL expansions, it follows that renewable generation /?'(£) = 
/(exp(Wj,(f,£)) is a function on the same stochastic space. 
We neglect the (small) noise in the conversion of wind speed 
into power and approximate / as a cubic spline interpolation 
through the filtered rated power data. This approximation is 
depicted by the red line in Fig. 1. 



Fig. 7. Decay of covariance matrix Ly components, shown with thin grey 
lines, with increasing time lag At = \tj — r,|, anchored at several time instances 
i for site #16238. Values are normalized by the diagonal entries. Matern 
covariance model is shown with blue line. 

TABLE II 

Parameters of Matern kernels corresponding to select wind 

sites. 


Wind Site 

It 

V 

15414 

11.40 

0.56 

16238 

11.15 

0.57 

3560 

9.79 

0.78 


B. Wind Power Forecast Models 


In this section we discuss the algorithm for generating wind 
power scenarios for day-ahead economic dispatch studies. 
We employ the KLE approach described in the previous 
section, with the covariance matrix adjusted to account for 
typical uncertainties of wind speed data obtained from weather 
forecast moedls. 

We first start by fitting a functional form through the 
covariance matrix corresponding to each site. Fig. 7 shows, 
with thin grey lines, a typical decay in the components of Eyy 
with increasing time lag. We model this data with a Matern 
covariance kernel [18] 


- /A , 2 1-v ( V2vAr\ ( V2vA t\ 

M—) 


( 8 ) 


Here, v and l t are positive parameters, T(-) is the gamma 
function, and K v (-) is the modified Bessel function of the 
second kind. The Matern kernel offers more flexibility for 
modeling covariances compared to, for example, exponential 
and square exponential forms which are particular cases of the 
Matern kernel, for v = 1/2 and v —► 00 , respectively. Table II 
shows the parameters of the Matern kernels for the three 
sites employed in this study. We obtain similar values for the 
two sites (15414 and 16238) that are geographically close. 
Furthermore, the values for the v parameter for these two sites 
are close to v = 1/2 indicating covariance matrices that are 
close to exponential form. Next we estimate the magnitude of 
the variance in the estimates of wind magnitude. Since data 
on wind forecast and actual realizations were not immedi- 
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ately available to us, we proceeded to estimate uncertainty 
in the wind forecasts given available data for wind power. 
Specifically, we used data for day-ahead forecast and actual 
realizations obtained from Belgium Electricity Grid Operator 
ELIA [19]. We considered the aggregate wind power values 
for the land-based wind farms and computed the standard 
deviation for hourly wind power output based on data for 
years 2012 through 2015. We found a relative value of about 
Op = 35%. Further, this value is independent of the time of 
the day for the day-ahead forecasts. 

Based on the formulations presented in this section we em¬ 
ploy the following algorithm to generate wind power samples 
that are consistent to historical wind characteristics at the 
selected wind sites: 

1) Forecast day-ahead wind profiles at selected wind sites. 

2) Using Op defined above, estimate a mean Ow based 
on the mean forecast wind speed and the rated power 
output curves for each site, similar to the one presented 
in Fig. 1. 

3) Construct covariance matrix, E w = C^Ejy, where E^ is 
the normalized Matern covariance expression presented 
in Eq. (8). 

4) Perform eigen-decomposition for Ew and generate wind 
samples via Eq. (5) for select samples of 

5) Convert wind samples into wind power values, see 
Fig. 1. 

IV. Accurate estimation of expected cost with 

LIMITED SAMPLES 

We now return to the evaluation of the expected cost 
in Eq. (1) corresponding to the ED problem in Eqs. (2) 
and (3). We can estimate the expected production cost by 
using a finite number of renewable power realizations (i.e., 
scenarios) s € § sampled from the joint density PDF( £). For 
our current example, each ^ is a standard normal RV, hence 
the sampling can be done independently in each stochastic 
direction. Defining p = 1/|S|, where |S| is the cardinality of 
§, the stochastic ED in Eq. (1) can be rewritten as: 

min p£Q(*A) (9) 

/,p,9,0 iG§ 

where Ql'x.s) is the solution of Eqs. (2) and (3) for a particular 
instance of the renewable generation p r r ( i;) —> // r (s). 

Formulation (9) represents an extensive form of the stochas¬ 
tic ED problem, based on |S| sampled scenarios from the 
stochastic space corresponding to wind power generators. 
The typical scenario sampling approach described above uses 
Monte Carlo (MC) sampling to approximate an integration, 
thereby estimating an expectation. While MC algorithms are 
commonly used for their convenience and robustness, their 
poor convergence rate is well-known. The MC estimate of the 
expectation has error 

W[Q(xX)]/VWl (10) 

where V[Q] denotes the variance of the RV Q. Given the signif¬ 
icant additional complexity incurred by including stochasticity 
in the optimization problem, a stochastic formulation becomes 


advantageous relative to a deterministic formulation when the 
variance is large. Hence, accurate estimation of the expectation 
is not only an academic exercise but is important in practice. 

According to Eq. (10), accurate estimation can be achieved 
by increasing the number of samples. However, a linear 
decrease in error requires a quadratic increase in the number of 
samples, which can quickly render the stochastic optimization 
problem intractable. This illustrates the limitation of MC 
algorithms in providing accurate estimations; while they are 
convenient, they are not efficient. 

We propose a method based on Polynomial Chaos expan¬ 
sions that can enable high precision quantification of uncer¬ 
tainties with fewer samples. In the following sections, we will 
first outline the Polynomial Chaos expansion construction and 
then present its implementation for the ED problem. 


A. Representation of uncertainty using Polynomial Chaos 

Given the formulation in Eq. (3) with uncertain/random 
loads leading to uncertain/random production costs, we em¬ 
ploy efficient UQ methods that rely on functional represen¬ 
tations of random variables. Specifically, we use Polynomial 
Chaos (PC) expansions. A brief description of PC is presented 
below. For an in-depth description, the reader is referred to a 
series of publications on this topic [11], [20]-[22]. 

We begin by setting up a requisite theoretical framework 
as follows. Define the probability space (£2,©,P), where Cl 
is a sample space, © is a o-algebra on Cl, and P is a 
probability measure on (£2,6). Further, defining the germ 
\ = {^ 1 ,^ 2 as a set of independent identically dis¬ 
tributed ( iid ) RVs in L ?(£2, &,P), to be further specified below, 
we focus on the probability space (Cl,6^,P) employing the 
sigma algebra generated by In this framework, any RV 
X : Cl —► R, where by construction X £ L 2 (Cl,&^,P), can be 
written as a PC expansion (PCE): 


oo 

X(co) =X(ij(a>)) = £a* v P*(£) 


k =o 


( 11 ) 


where the basis functions 'P< are multivariate polynomials* 
that are orthogonal, by construction, with respect to the density 
of £. Thus 

W ’j) = I j(£,)dP(£,) = M'p?) (12) 


where 8, ; is Kronecker’s delta. Further, given this orthogonal¬ 


ity, we have 


oh = 


(ggt) 


(13) 


where the inner product is defined, for any RV Z(§), by the 
Galerkin projection 


(Z) = J (14) 

Moreover, the 'Pjt are products of univariate polynomials, 
namely ¥*(£) = t|4, (qi) • ■ ■ (^«)- In a practical computa¬ 

tional context, one truncates the PCE to order p. The number 


* Generally, other, non-polynomial basis functions can be used, but here we 
restrict ourselves, without loss of generality, to the most common polynomial- 
based usage. 
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of terms in the resulting finite PCE 

p 

& k {%) (15) 

k—0 

is given by P+1 = (n + p)\/n\p\. We dispense with the sa 
symbol in the remainder of this paper, employing for any RV 
X (4) its truncated PCE 

p 

X='£a k 'V k (Z). (16) 

k =0 

Generalized PC (gPC) expansions have been developed by 
Xiu and Karniadakis [22] using a broad class of orthogonal 
polynomials in the “Askey” scheme [23], Each family of 
polynomials corresponds to a given choice of distribution for 
the and is, by construction, orthogonal with respect to the 
density of In general, the most useful choices for (£, V P) 
are normal RVs with Hermite polynomials and uniform RVs 
with Legendre polynomials. 

B. Construction of PCE for the Minimum Cost 

In the context of this study with uncertain renewable wind 
generation dependent on a stochastic space of iid standard 
normal random variables, £ = • • • ,^ n }, we employ 

Hermite polynomials to construct a Hermite-Gauss (HG) PCE 
for the (minimum) cost Q(x. £). We employ Eq. (16) for a 
truncated HG PCE 

p 

Qpc{x,$) = 'Ec k (xye k (%), (17) 

k =0 

where 'P*(£) are n-variate Hermite polynomials. The coeffi¬ 
cients c k depend on the discrete variable x, hence separate PCE 
approximations for Q will be constructed for each instance of 
x chosen in the unit commitment stage. Given Eq. (14), we 
have 

c k (x) = ( fj* } = . ' [ Q{x£)'¥ k {%)p{£ > )dk. (18) 

where we have used p(ff) as a product of univariate standard 
normal PDFs. The dimensionality of the PCE in Eq. (17) is the 
same as the number of stochastic dimensions used to represent 
the uncertain renewable power. For the study presented in this 
paper, n = 16. 

Given <2 pc(*i £)> then, we have 

Q(x) = E, % \Q{x&)\ = (Q(x,\)) = co, (19) 

being the solution of the stochastic ED problem. The PCE in 
Eq. (17) can also be used to generated higher-order moments 
for the minimum cost, based directly on the coefficients for 
the corresponding PC basis terms. 

Several methods can be employed to evaluate the projection 
integrals in Eq. (18). MC methods can be used in principle, 
but are impractical given their slow convergence rate. Alterna¬ 
tively, for smooth integrands, and particularly in low-moderate 
dimensional problems, sparse quadrature methods [24]-[26] 
can provide highly accurate results with smaller numbers of 
deterministic samples. Fig. 8 shows, in a 2D configuration, the 
locations of deterministic samples we use, with a sparse grid 
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Fig. 8. Placement of deterministic samples via a sparse grid approach: (a) 
Level 1, (b) Level 2, (c) Level 3, and (d) Level 4. For Levels 2 through 4, 
the red symbols show samples for previous levels while blue symbols show 
newly added samples. 


employing Gauss-Rronrod quadrature [27]. Several levels are 
shown in the figure, starting with the first level in Fig. 8(a), 
and following with additional samples leading to Levels 2, 3, 
and 4, in Fig. 8(b)-(d). A first order PCE requires a Level 
2 sparse grid, while a second order PCE requires a Level 3 
grid. The number of requisite samples using sparse-quadrature 
evaluation of the projection integrals, for a given requisite 
surrogate accuracy, is much smaller than the corresponding 
number of MC samples, as we will illustrate in the next 
subsection. 

C. Numerical Results 

We now present numerical results comparing the expected 
costs of our SED model using scenarios obtained with both our 
Polynomial Chaos approach to sampling versus a traditional 
Monte Carlo approach. We consider the IEEE 118-bus test 
system [28] augmented with the three wind generation sites 
discussed in Section III. The three renewable generators at 
sites 15414, 16238, and 3560 replace the conventional gener¬ 
ators at buses 89, 69, and 10, respectively. We employ |7j = 24 
time periods and represent stochastic renewable generation via 
the Karhunen-Loeve expansions presented in Section III. We 
use 6 KL modes per site. Given the strong dependence between 
the first two RVs for sites 15414 and 16238, the effective 
dimensionality of the system decreases by 2, to a total of 16. 

We first explore the dependence of Q(x,^) on the stochastic 
space that characterizes renewable generation. Fig. 9 shows 
2D slices through the 16D \ space. All slices are anchored 
at the origin £ = 0. RVs qi through qo correspond to site 
16238, while qi, and ^7 - 4io correspond to site 15414. 
The remaining random variables, <qn - q !6 , correspond to 
site 3560. The expected cost is about $2.4 million, and the 

















Fig. 9. Two dimensional slices through the center of the stochastic space, 
£ = 0. Random variables corresponding to (a) first two modes for Wyoming 
sites, (b) 1st and 3rd mode for WY sites, (c) first modes for WY and CA 
sites, and (d) first two modes for CA site. The countour lines correspond to 
iso-levels of Q, increasing from blue to yellow to red. 


relative change about the mean is about 15% in the numerical 
tests presented here. The results shown in the top row of 
Fig. 9 indicate that Q is strongly dependent on ^i, while 
the dependence on 'tp_ and c,j is weaker - as suggested by 
negligible contour changes in those directions. The slice in 
the (£,i, £,ii) plane highlights the contribution of to the 
variation of Q, while the slice in the (^ 11 ,^ 12 ) plane indicates 
less impact of the second mode for site 3560 on cost. Other 
2D slices (not shown) confirm the diminishing impact that 
higher-order modes have on Q. Collectively, these results also 
indicate that Qix.E,) is smooth in £,, which makes it amenable 
to a PCE representation. 

We next proceed to test the accuracy of the PCE represen¬ 
tation for Q(x, 4) with respect to actual model evaluations. 
We construct several PCEs, from 1st to 4th order. Given the 
16-dimensional stochastic space corresponding to the three 
wind sites, the number of sparse quadrature sample points 
necessary for the construction of the PCE coefficients is 
approximately 5 x 10 2 for a 1st order expansion, 5 x 10 3 for a 
2nd order expansion, 3.6 x 10 4 for a 3rd order expansion, and 
2 x 10 5 for a 4th order expansion. We cross-validate our PCE 
representations relative to exact solutions at 5 x 10 5 randomly 
chosen £ samples. Fig. 10 shows the relative L\ error between 
actual model evaluations and the PCE-approximated cost. This 
error is relative to the expected minimum cost computed using 
all model evaluations available and is converted to percentages. 
For brevity, only a random subset of 2 x 10 4 samples are 
shown. These results indicate a drop of about one order of 
magnitude, per polynomial degree, for the magnitude of the 
error between full model evaluations and the PCE results. 

The above analysis explores the dependence of minimum 



sample # 

Fig. 10. Cross-validation of PCE approximation for the minimum cost 
Q(x, lj). The blue, red, green, and black symbols correspond to 1st through 
4th order PCEs. 


cost for individual samples on the stochastic space correspond¬ 
ing to the uncertain wind at select sites, and the accuracy 
of the PCE approach in capturing this dependence. However, 
in the context of the SED problem, and for the present 
work in particular, we are only interested the accuracy of the 
expected minimum cost (Q(x .^)).' : Hence, our goal here is to 
demonstrate the efficiency of a quadrature approach to estimate 
the expectation. 

Toward this goal, we show the convergence of (Q{x,%)) as 
a function of the number of samples considered. For the PCE 
approach, the expected minimum cost is simply given by the co 
coefficient in Eq. (18). The results shown with blue are based 
on this coefficient, computed with increasingly accurate sparse 
quadratures. The error measure is relative, being normalized 
with respect to the value corresponding to the next higher 
accuracy result, defined as 


Epc,i 


ko,--co,- +1 | 

Gh'+i 


( 20 ) 


where the subscript i denotes the quadrature level, e.g., 33 
model evaluations for level 1,513 model evaluations for level 
2, and so forth. The error for the MC results is computed 
using a similar approach. However, because results for the 
MC approach depend on a randomly drawn set of samples, 
we show results for several realizations. The specific formula 
is given as 


F J — 
c MCA ~ 


\Qi-Q l+ A 

Qi +1 


( 21 ) 


where the subscript i denotes the number of samples used to 
compute the average cost Q t , e.g., 10 samples for i = 1, 100 
samples for i = 2, etc. We employed 10 realizations, indexed 
by superscript j, for each set of MC samples; the average 
cost over all realizations is denoted using a double overline. 
Eha is shown with red circles in Fig. 11, while the red line 
connects the mean error over all realizations for a particular 
number of MC samples. The dashed lines show power-law 
curve fits, a x N~ b , for the data, which we use to analyze 
convergence rates as a function of the number of samples. 
For the MC approach we recover the theoretical convergence 
rate, b = 0.5, while for the PCE approach b = 1.9. Visual 


1 Our approach is not restricted to expected costs; other moments of Q can 
be similarly considered. 
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inspection of the PCE results and power law fits indicate a 
stronger polynomial convergence rate than MC, rather than an 
exponential rate (results not shown) corresponding to spectral 
accuracy. 

In this figure a 1% error level corresponds to about $24K, 
while a 0.01% error level is in the hundreds of dollars. These 
results indicate that if a low degree of accuracy, e.g. 1%, 
is sufficient or desired, then the PCE approach exhibits a 
computational cost that is comparable to the MC approach. 
Nevertheless there is a considerable spread in the accuracy of 
these estimates, making MC estimates with a small number 
of samples somewhat unreliable. For situations where higher 
accuracies are required, for example of 0.01%, the PCE 
estimates require one to two orders of magnitude less model 
evaluations compared to the MC estimates. 



Fig. 11. Comparison of errors in the estimation of expected cost between 
PCE-based results and Monte Carlo estimates. Dashed lines show power-law 
fits, axN b to test the convergence rates of the MC and PCE approaches. 


V. Conclusion 

In this paper, we present methods for efficient representation 
of uncertainty, with emphasis on the Stochastic Economic Dis¬ 
patch problem. We present two, related, category of methods, 
one for spectral random fields representations, and one for 
functional representation of random variables. 

In the first part of this paper we employ data on renewable 
wind generations provided by NREL. We determine that 24- 
hour wind samples can be represented via Karhunen-Loeve 
(KL) expansions. This expansion represents stochastic pro¬ 
cesses by a linear combination of orthogonal modes. Further, 
the KL representation employs the basis that provides the 
minimum total L 2 error. Analysis of the KL eigenspectra at 
several wind sites indicate that 6 KL modes are sufficient to 
capture about 95% of the total variance, effectively reducing 
the dimensionality of the stochastic space by a factor of 4. 
The dimensionality of the stochastic space is further reduced 
by exploiting the dependence between KL random variables 
for wind sites that are geographically close. 

In the second part of the paper we present an approach to 
reduce the computational cost associated with stochastic unit 
commitment and economic dispatch, by reducing the number 


of required forecast samples. This approach is based on 
Polynomial Chaos expansion (PCE) models for the production 
cost that cover the uncertainty in the renewable generation. The 
construction of the PCE terms is based on the projection of 
the model on increasingly higher basis modes. Consequently, 
the global error in an L 2 sense between the surrogate model 
and the actual simulations is easily controlled. 

We present computational results for the 118-bus test case, 
augmented to account for renewable generation. We find 
that quadratic PCE models for the production cost showed 
pointwise relative errors less than 1% throughout the uncer¬ 
tain demand space. For relative accuracies of 0 (10 4 ), the 
PCE approach is one to two orders of computationally more 
efficient compared to Monte Carlo estimates. 

We plan to extend the methods presented here to higher¬ 
dimensional power grid problems. In order to alleviate the 
curse of dimensionality in such studies we plan to emply 
global sensitivity analysis to eliminate the stochastic variables 
that are not important for the expected cost. And based on 
results observed in this paper, we plan to explore adaptive 
sparse quadrature constructions to tailor of the PCE order to 
the specific dependence of minimum cost on each component 
of the stochastic space. 
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Appendix. Total Variance of Karhunen-Loeve 
Expansion 

Starting from Eq. (5), it follows that 

' 24 

Var[W L (f,C 0 )] = Var £ s/X k f k {t)t, k 

_k—\ 

Here, we converted the infinite sum into one over the 24 modes 
that can be extracted from 24h wind samples. As noted in 


Section III, ^ are uncorrelated random variables with zero 
mean and unit variance. To simplify notation, in this Appendix, 
we change notation and denote the product \fX k ^ k as A/,. This 
new set of RVs are uncorrelated with zero mean and variance 

X k . 

The variance of the rhs of the expression above can be 
written as 


Var Y,Mt)Xk 
k 


= E £/t(f)X*-E £/*(*)** =E £/*(0(XI-EN) 

V k k ) V k 


= E Y J fht)(X k -E[X k })- + 2Y J Mt)f j (t)(Xi-nXi})(Xj-E[Xj}) 

. k ¥j 

= I fk(t) E [(X k - E [A,]) 2 ] +2 ^Mt)fj(t) E [(A, - E [A,]) (Xj - E [A,-])] 


For uncorrelated random variables A,, Ay, the expecta¬ 
tion in the second sum in Eq. (24) is identically zero, 
E [(A; — E [A,-]) (Ay — E [Ay])] = 0. Making use of Var[A A ] = X k , 
it follows that 

Var £/*(r)*t =£/^(0Var[A A ]=£/ 2 (0^ (25) 

k k k 

The total variance over the entire deterministic space is given 




Given that modes f k are orthonormal over the deterministic 
space, total variance of the KL expansion is given by 



[ W L (t,a>)\dt 


Y*Xk 



